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Abstract 

TSIL is a library of utilities for the numerical calculation of dimcnsionally regularized two-loop 
self-energy integrals. A convenient basis for these functions is given by the integrals obtained at 
the end of O.V. Tarasov's recurrence relation algorithm. The program computes the values of all 
of these basis functions, for arbitrary input masses and external momentum. When analytical 
expressions in terms of polylogarithms are available, they are used. Otherwise, the evaluation 
proceeds by a Runge-Kutta integration of the coupled first-order differential equations for the 
basis integrals, using the external momentum invariant as the independent variable. The starting 
point of the integration is provided by known analytic expressions at (or near) zero external 
momentum. The code is written in C, and may be linked from C/C++ or Fortran. A Fortran 
interface is provided. We describe the structure and usage of the program, and provide a simple 
example application. We also compute two new cases analytically, and compare all of our 
notations and conventions for the two-loop self-energy integrals to those used by several other 
groups. 
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1 Program summary 

Title of program: TSIL 
Version number: 1.1 



Available at: http://zippy.physics.niu.edu/TSIL 



http : //www . otterbein . edu/home/f ac/dvdgrbrt/TSIL 
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Programming Language: C 

Platform: Any platform supporting the GNU Compiler Collection (gcc), the Intel C compiler (ice), 
or a similar C compiler with support for complex mathematics 

Distribution format: ASCII source code 

Keywords: quantum field theory, Feynman integrals, two-loop integrals, self-energy corrections, 
dimensional regularization 

Nature of physical problem: Numerical evaluation of dimensionally regulated Feynman integrals 
needed in two-loop self-energy calculations in relativistic quantum field theory in four dimensions 
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Method of solution: Analytical evaluation in terms of polylogarithms when possible, otherwise 
through Runge-Kutta solution of differential equations 



Limitations: Loss of accuracy in some unnatural threshold cases that do not have vanishing masses 
Typical running time: Less than a second 

2 Introduction 

The precision of measurements within the Standard Model requires a level of accuracy obtained 
from two-loop, and even higher-order, calculations in quantum field theory. The future discoveries 
of the Large Hadron Collider and a future e + e~ linear collider may well extend these requirements 
even further. For example, supersymmetry predicts the existence of many new particle states with 
perturbative interactions. The most important observables in softly broken supersymmetry are just 
the new particle masses. Therefore, comparisons of particular models of supersymmetry breaking 
with experiment will require systematic methods for two-loop computations of physical pole masses 
in terms of the underlying Lagrangian parameters. 

In this paper, we describe a software package called TSIL (Two- loop Self-energy Integral Li- 
brary^ that can be used to compute two-loop self-energy functions. In a general quantum field 
theory, the necessary dimensionally regularized Feynman integrals can be expressed in the form: 

„8-2d f f (fc 2 r( g 2 r(fc-pr3( g .p)"4 (fc . g) n 5 

^ J J q [k 2 +x]^[q 2 + y]^[{k-pY + z]^[{q-p) 2 +u]^[{k-q) 2 + vY^ 1 ' ' 

for non-negative integers m, r{. (See the master topology M in Figure[TJ) The integrations are over 
Euclidean momenta in 

d = 4-2e (2.2) 

dimensions. An algorithm has been derived by O.V. Tarasov [I], and implemented in a Mathematica 
program TARCER by Mertig and Scharf [2], that allows all such integrals to be systematically 
reduced to linear combinations of a few basis integrals. The coefficients are ratios of polynomials in 
the external momentum invariant and the propagator squared masses x, y, z, u, v. In the remainder 
of this section, we will describe our notations and conventions |] for the two- loop basis integrals and 
some related functions, and describe the strategy used by TSIL to compute them. 
First, we define a loop factor 

C =^ 2 )^y = (^A 2 . (2.3) 

The regularization scale [i is related to the renormalization scale Q (in the modified minimal sub- 
traction renormalization scheme based on dimensional regularization [6], or dimensional reduction 



'In the Hopi culture indigenous to the American southwest, Tsil is the Chili Pepper Kachina, one of many 
supernatural spirits represented by masked doll-like figurines and impersonated by ceremonial dancers. Tsil is one of 
the runner Kachinas. When he overtakes you in a race, he may stuff your mouth with hot chili peppers. 

'These are the same as in refs. [4j[5]; the correspondences with some other papers is given in Appendix A. 
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Figure 1: Feynman diagram topologies for the one- and two- loop vacuum and self-energy integrals 
as defined in this paper. 



[7] for softly-broken super symmetric theories) by 

Q 2 = 4vre-V. 

Logarithms of dimensionful quantities are always given in terms of 

liLY = ln(X/Q 2 ). 

The loop integrals are functions of a common external momentum invariant 

s = —p 2 



(2.4) 



(2.5) 



(2.6) 



using a Euclidean or signature (—+++) metric. (Note that the sign convention is such that for 
a stable physical particle with mass m, there is a pole at s = m 2 .) On the physical sheet, s has 
an infinitesimal positive imaginary part. Since all functions in any given equation typically have 
the same s and Q 2 , they are not included explicitly in the list of arguments in written formulas. 
A prime on a squared-mass argument of a function indicates differentiation with respect to that 
argument. Thus, for example 

" d 2 



f(x',x,z') 



dydz 



f(y,x,z) 



J y=x 



We now define one-loop vacuum and self-energy integrals as: 

1 



A(x) 



C / d a k 



[k 2 + x] 



1 



B(x ' y) = c I ddk WTx][(k- P y + y] 

We also define two-loop integrals according to: 

S(W) = ° 2 1 ^ I ^ W 2 + A [Q 2 +V][(k + q- P) 2 + z] 



(2.7) 

(2.8) 
(2.9) 

(2.10) 
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I(%,y,z) = S(x,y,z)\ s=0 
T(x,y,z) = -S{x',y,z) 
V(x,y,z,u) = C 2 J d d k J d d q 



[k 2 + x][(k - p) 2 + y][q 2 + z][(q + k- p) 2 + u] 



V(x,y,z,u) 
M(x,y,z,u,v) 



-U(x,y',z,u) 
C 2 J d d k J d d q- 



(2.11) 
(2.12) 
(2.13) 
(2.14) 

r- ( 2 -i5) 



[k 2 + x] [q 2 + y] \{k - p) 2 + z] [(q - p) 2 + u] [(k - q) 2 + v] 

The corresponding Feynman diagram topologies are shown in fig. [TJ These integrals have various 
symmetries that are clear from the diagrams: B(x,y) is invariant under interchange of x,y; the 
"sunrise" integral S(x,y, z) and I(x,y, z) are invariant under interchange of any two of x,y,z; 
T(x,y,z) is invariant under y <-> z; XJ(x,y,z, u) and V(ar, y, z, u) are invariant under z <-> u; and 
the "master" integral M(x, y, z, u, v ) is invariant under each of the interchanges (x,z) <-> (y,u), 
and (x,y) <-> (z,u), and «-> (u,z). 

It is often convenient to introduce modified integrals in which appropriate divergent parts 
have been subtracted and the regulator removed. At one-loop order, we define the finite and 
e-independent integrals: 



A(x) 
B(x,y) 
At two loops, we let 

S(x,y,z 

where 



lim [A (a;) + x/e] = x(lnx — 1) 



lim [B(x,y) - 1/e] 

e— >U 



dt ln[te + (1 - t)y - t(l - t)s). 



lim 



S(x, y, z) - S£(ar, y, z) - S£l(x, y, z) 



(2)/ 



(A(x)+A(y) + A(z))/e 

(x + y + z)/2e 2 + (a/2 -x-y- z)/2e 



(2.16) 
(2.17) 

(2.18) 



(2.19) 
(2.20) 



are the contributions from one-loop subdivergences and from the remaining two-loop divergences, 
respectively. In addition, 



I(x,y,z) 
T(x,y,z) 



Similarly, we define 



where 



U(x,y,z,u) 



lim 



S(x,y,z)\ s=0 
-S{x',y,z). 

(1), 



U(x, y ,z,u)-U£>(x,y)-U£l 



U {2) 



B(x,y)/e 
-l/2e 2 + l/2e 



(2.21) 
(2.22) 



(2.23) 



(2.24) 
(2.25) 
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are again the contributions from one-loop sub-divergences and the remaining two-loop divergences. 
Also, we define 

V(x,y,z,u) = -U(x,y',z,u). (2.26) 
The master integral is free of divergences, so we define 

M(x,y, z,u,v) = limM(x, y, z,u,v). (2.27) 

e— >0 

Thus, bold-faced letters A, B, I, S, T, U, V represent the original regularized integrals that di- 
verge as e — > 0, while the ordinary letters A, B, I, S, T, U, V, M are used to represent functions that 
are finite and independent of e by definition. Note, however, that as defined above /, S, T, U, V are 
not simply the e-independent terms in expansions in small e. The following expansions are useful 
for converting between I, S, T, U, V and I, S,T,U,V: 

A(x) = -x/e + A(x) + eA e (x) + 0(e 2 ) (2.28) 

B{x,y) = l/e + B(x,y) + eB e (x,y) + 0{e 2 ), (2.29) 

where 

A t (x) = x[-l-C(2)/2 + Ex- (Imc) 2 /2] (2.30) 

B e {x,y) = C(2)/2 + ^ f 1 dt(Htx + {l-t)y-t{l-t)s}) 2 . (2.31) 

2 Jo 

Here £ is the Riemann zeta function. The function B e (x, y) can be expressed in terms of dilogarithms 
[8], and is given by the coefficient of 5 in eq. (83) of ref. [9]. Also, 

B(x', y) = [(3 - d)(s - x + y)B{x, y) + (2 - d){A{y) + (s - x - y)A(x)/2x}] /A sxy (2.32) 

where 

A abc = a 2 + b 2 + c 2 - 2ab - 2ac - 26c. (2.33) 

From the preceding equations it follows that 

I(x,y,z) = -{x + y + z)/2e 2 + [A{x)+A{y)+A{z)-{x + y + z)/2}/e (2.34) 

+I(x, y, z) + A e (x) + A e (y) + A e {z) + O(e) 

S(x,y,z) = -(x + y + z)/2e 2 + [A(x)+A(y)+A(z)-(x + y + z)/2 + s/A}/e (2.35) 

+S(x, y, z) + A e (x) + A e (y) + A e {z) + O(e) 

T(x,y,z) = l/2e 2 - [A(x)/x + 1/2] /e + T(x,y,z) + [A(x) - A e (x)]/x + 0(e) (2.36) 

U(x,y,z,u) = l/2e 2 + [1/2 + B(x,y)] /e + U(x, y,z,u)+B e (x, y) + 0{e) (2.37) 

V(x,y,z,u) = -[(s + x-y)(B{x,y)-l) + 2A(x) + (s-x-y)A(y)/y]/A sxy (2.38) 



+V(x, y, z, u) + {(s + x - y)[B e (x, y) - 2B{x, y)] + 2A e (x) - 2A{x) 
x - y)[A e (y) - A(y)]/y}/A sxy + 0(e) 



+ s 



M(x,y,z,u,v) = M(x,y,z,u,v) +0(e). (2.39) 



G 



The internal workings of the TSIL code use the functions A, B, I, S,T,U,V, M rather than 
their bold-faced counterparts. This is because we find that renormalized expressions for physical 
quantities are more compactly written in terms of the non-bold-faced integrals. However, both 
types of functions are available as outputs, with the proviso that for I, S, T, U, V, M we keep only 
the pole and finite terms as e — > 0, and for A, B only terms through order e. 

Tarasov's algorithm [JJ allows any integral of the form of eq. (|2.ip to be expressecjl as a linear 
combination of the following two-loop basis integrals: 

M(x,y,z,u,v), U(z,x,y,v), U(u,y,x,v), U(x,z,u,v), U(y,u,z,v), T(v,y,z), 
T(u,x,v), T(y,z,v), T(x,u,v), T(z,y,v), T(v,x,u), S(v,y,z), S(u,x,v), (2.40) 

plus terms involving the two-loop vacuum integrals I(x,y,v) or I(z,u,v), or quadratic in the one- 
loop integrals. 

In particular, the V and V integrals can be expressed as linear combinations of the others (see 
eqs. (3.22)^(3.28) and (6.21) of ref. [1]), and so are not included in the basis. However, they are 
included as outputs, because some results are more compactly written in terms of them. Because 
T(x, y, z) is divergent in the limit x — > 0, it is also sometimes useful to define the function: 

T(x, y, z) = T(x, y, z) + B(y, z) Ex. (2.41) 

For x = 0, this function is well-defined and can be written in terms of dilogarithms (see eqs. (6.18)- 
(6.19) of ref. [4]). In that limit it can also be rewritten in terms of the other basis functions, see 
eqs. (A.15)-(A.16) of ref. [5], but is still available as an output of TSIL for convenience. 

It remains to provide a means for the numerical computation of the basis integrals. For special 
values of masses and external momentum, it is possible to compute the two-loop integrals analyt- 
ically in terms of polylogarithms |10| . This requires |llj that there is no three-particle cut of the 
diagram for which the cut propagator masses m\, 777-2, w*3 and the five quantities 

Scut, Scut - (mi ± 777,2 ± m 3 ) 2 , (2.42) 

(where s cu t = — p^ nt is the momentum invariant for the total momentum flowing across the cut) 
are all non-zero. Many analytical results for various such special cases have been worked out in 
refs. |12j-|20j. [9], [I], and Appendix B of the present paper. There are also expansions |21|- 
|25j in large and small values of the external momentum invariant, and near the thresholds and 
pseudo-thresholds [26]-[35]. Analytical results in terms of polylogarithms for the A, B, I, S, T, U, M 
functions at generic values of s are reviewed in section VI of ref. [I] . These and some other analytical 
formulas for special cases have been implemented in the TSIL code. They consist of the master 

§ Actually, Tarasov's algorithm relates the general integral to the bold-faced versions of the basis functions, and 
holds for general d. By neglecting contributions that vanish for e —* 0, one can write the results in terms of the 
non-bold-faced functions. 
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integral cases: 



M(x,x,y,y,0), M(0, 0, 0, x, 0), M(0,x,0,x,x), ref. [H] 

M(0,0,0,0,x), M(0, x, 0, y, 0), ref. [9] 

M(0,0,0,x,x), M(0,2J,a:,0,0), M(0,x,x,x,0), ref. pi] 

M(0, 0, x, y, 0), Appendix B 

M(0,2J,a:,0,a;)| s=a; , ref. pi] 

M(0,y,y,0,x)| s=x , Appendix B 

M(0,x,y,0,y)\ s=x , refs. [36], [37] (see also Appendix B) 

and functions obtained by permutations of the arguments, and all subordinate integrals A, B, I, 

S, T, U, V, T that have topologies obtained by removing one or more propagator lines from the 
cases above. These include: 



U(x,y,0,y), ref. [E] 

U(x,y,0,0), [7(0,0, 0,x), ref. [9] 

U(0,x,y,z), ref. 0] 

U(x,0,0,y), ref. [5] 

U(x, 0, y, y)\ s=x , U(y,0,y,x)\ s=x , refs. [38] [39] (see also Appendix B) 

and 

S(0,x,y), T(x,0,y), T(0,x,y), ref. [IS] 

S'(^,2/,2/)|s=a : ) r(z,y,y)| s=a; , T(y, x, y)| s=x . refs. [38J [29] 

Also included are all of the functions at s = 0, which can be easily expressed in terms of the A 
and / functions and derivatives of them, which are in turn known [141 |4T)] analytically in terms of 
logarithms and dilogarithms. 

For the case of generic masses and s, another method is needed. Integral representations have 
been studied in refs. [H]-[50]. For TSIL, we instead use the differential equations method [5T]-[52] 
to evaluate the integrals numerically. For the S, T and U integrals, this strategy was proposed and 
implemented in [S3] _ [SZ]- The method was rewritten in terms of the S, T, U integrals and extended 
to M in ref. [3]. To see how the method works, consider the functions listed in eq. (|2.40p and also 
B{x,z) and B(y,u) and the product B(x, z)B(y,u). Let us denote these sixteen quantities by F{ 
where i = 1, . . . , 16. Considered as functions of s for fixed x, y, z, u, v, Q 2 they can be shown to 
satisfy a set of coupled first-order differential equations of the form 

^- s Fi = J2 C iJ F j+ C ^ ( 2 - 43 ) 
j 

where the coefficients Cij and Cj are ratios of polynomials in the squared masses and s (and of 
the analytically known A and / functions, for some of the coefficients not multiplying two- loop 
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Im[s] 



III | Re[s] 

thresholds 

Figure 2: When s is greater than or equal to the smallest non-zero threshold or pseudo-threshold, 
the Runge-Kutta integration proceeds along a contour with the shape shown here, as suggested in 
ref. [56j. 

functions). The coefficients can be evaluated by using Tarasov's recurrence relations, and were 
listed in [3J. 

At s = the values of all of the functions and their derivatives with respect to s (and/or their 
expansions in small s) are known analytically in terms of dilogarithms. Therefore, one can integrate 
the functions^) by the Runge-Kutta method to any other value of s. In order to avoid numerical 
problems from integrating through thresholds and pseudo-thresholds, we use the suggestion of 
ref. [56] by following a displaced contour in the upper-half complex s plane, as shown in fig. [21 
whenever s is greater than the smallest non-zero threshold or pseudo-threshold. This contour 
starts from s = (or, in some cases, a point very close by, as explained in the next section) and 
proceeds to a point is- im , from there to s + zsjrn, and from there to the desired value s. Here Sim is 
real and positive. Since s has an infinitesimal positive imaginary part on the physical sheet, this 
procedure also automatically produces the correct branch cut behavior for functions when s is above 
thresholds. When s is less than or equal to the smallest non-zero threshold or pseudo-threshold, 
the Runge-Kutta integration instead proceeds directly along the real s axis. In typical physical 
problems, if the master integral is needed, then so will be all of its subordinate B, S, T, U integrals. 
These are all obtained simultaneously as a result of the Runge-Kutta method. Furthermore, checks 
on the numerical accuracy can be made by varying the Runge-Kutta step size parameters and the 
choice of contour in the upper-half complex s plane. 

Most of the practical difficulties in realizing this program have to do with numerical instabilities 
when the final value of s is at or near a threshold or pseudo-threshold, or when the starting point 
s = is itself a threshold. We describe these issues and the methods used by TSIL to successfully 
evade them in the next section. 

^For the master integral, we actually run sM(x,y, z,u,v). 
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3 Numerical considerations near thresholds and pseudo-thresholds 



The two-loop master integral function M(x, y, z, u, v), and its subordinates listed in eq. (|2,40p . have 
two-particle and three-particle thresholds at s equal to: 

t XZ = {VX + VZ) 2 , t yu = (y/y + y^) 2 ', (3.1) 

t X uv = {Vx + Vu + Vv) 2 , t yzv = (y/y + \[z + Vv) 2 - (3.2) 

At the thresholds, the integral functions have branch cuts, and they therefore develop imaginary 
part contributions for s > Sthresh- The pseudo-thresholds occur at s equal to: 

Pxz = (Vx - Vz) 2 , p yu = {^fy - Vu) 2 , (3.3) 

Pxuv = (~y/x + Vu + Vv) 2 , Puxv = {Vx- Vu + Vv) 2 , (3.4) 

Pvxu = (Vx + Vu- Vv) 2 , Pyzv = (-VV + Vz + Vv) 2 , (3.5) 

Pzyv = (VV ~ Vz+ Vv) 2 , Pvyz = (VV + Vz ~ Vv) 2 ■ (3.6) 

The integral functions are analytic at the pseudo-threshold points (unless they coincide with a 
threshold), but the coefficients in the differential equation have pole singularities at both the thresh- 
olds and pseudo-thresholds. For values of s close to both types of special points, one must therefore 
be careful in numerical evaluation to avoid undefined quantities and round-off errors. In this section, 
we discuss the ways in which TSIL avoids these problems. 

First, we consider the case that the initial point of the Runge-Kutta integration, s = 0, is 
actually a threshold. This occurs for master integrals M(0, y, z, 0, 0) and M(0, y, 0, u, v), and per- 
mutations thereof. In these cases, some of the basis integral functions have logarithmic singularities 
and/or branch cuts at s = 0. To deal with this, we make a change of independent variable to 

r = E(-s), (3.7) 

and start the integration at a point vq = ln(— iS) = ln(5) — iir/2, with 5 real, positive, and extremely 
small (of order the relative error of the computer arithmetic). The initial values of the integrals are 
obtained from the expansions in small s, given in section V of ref. [3]. This change of variables is 
also used when s = is close to, but not exactly a threshold (with the exact criteria adjustable by 
the user). This variable r is used for the first leg of the contour in the Runge-Kutta integration. 

Next, consider the case that the final value of s is at, or very near, a threshold s t h re sh- I n this 
case, we find that the change of variable 

r = ln(l - s/sthresh) (3.8) 

is effective, and so is used by TSIL for the final part of the Runge-Kutta integration. When the 
final value of s is exactly equal to a threshold, then the endpoint of the running is taken to be 
r = ln(<5) — ivr/2, with 5 again taken to be real, positive, and extremely small. 

We next describe the Runge-Kutta algorithms used, since they have some slightly unusual 
properties dictated by the need for control of precision near thresholds and pseudo-thresholds. 
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Consider a vector of quantities F that satisfy coupled first-order differential equations dF/dt = 
f(t, F). The general form for an explicit m-stage Runge-Kutta routine advancing the solution from 
t to t + h is given by: 

m 

F(t + h)=F(t) + hY / b i k i , (3.9) 
i=l 

where 

i-l 

k i = f(t + Cih,F(t)+h^2aijkj) (for i = 1, . . . ,m). (3.10) 
Here dy, 6j, and Cj are known as Butcher coefficients, and satisfy c\ = 0, aio = 0, and 

i-l 

Cj = otjj (for z = 2, . . . , m) (3-H) 
3=1 

in 

= 1 (3.12) 

i=l 

plus other, non-linear, constraints |58j . The algorithm is said to be of nth order if the error is 
proportional to h n+1 for sufficiently small h. In order to implement automatic step-size adjustment, 
we use a 6-stage embedded Runge-Kutta [59] with coefficients given by Cash and Karp [60]. These 
give not only a 5th-order step as in eq. (|3.9p . but a 4th-order step estimate of the same form with 
different coefficients b*. This gives an error estimate for each dependent variable, for each step: 

m 

AF(t + h) = h^2(bi - bt%. (3.13) 
i=i 

The theoretical estimate of the step size needed so that the error for each dependent variable is less 
than 5p times the increment of that variable is then given by: 

1/4 



S P 



(3.14) 



^Max(|AF| 
where S is a safety factor less than unity. 

However, in the present application there is a special problem because the final destination point 
might be equal to (or close to) a threshold or pseudo-threshold. There, the individual coefficients in 
the derivatives of the functions might be ill-defined (or subject to large numerical round-off errors, 
because of small denominators in individual terms), even though the basis functions themselves 
are well-defined. To avoid this problem, we instead need to use an m-stage Runge-Kutta with the 
crucial property q < 1 for all m, so that no derivatives are ever evaluated at the endpoint. There 
are no 4-stage, 4th-order solutions to the Butcher coefficient conditions with this property, but 
there are many 5-stage, 4th-order solutions. We chose, somewhat arbitrarily, the set: 

Cl = (0, 1/4, 3/8, 1/2, 5/8) (3.15) 

bi = (-1/15, 2/3, 4/3, -10/3, 12/5) (3.16) 

021 = 1/4, a 32 = 3/8, a 42 = 1/2, a 52 = 35/72, a 54 = 5/36, (3.17) 

dij = for other i,j. (3.18) 
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This procedure is not as efficient as the Cash-Karp 6-stage, 5th-order algorithm under normal 
circumstances, and does not provide an error estimate, so it is used only where needed for the very 
last Runge-Kutta step. 

The TSIL implementation of the coefficients Cij and Cj in eqs. (|2,43p is also done in a special 
way to avoid roundoff errors near thresholds and pseudo-thresholds. The expressions as given in 
[5] for these coefficients appear to have double (or higher) poles in s for some special values of 
the masses with degenerate thresholds and pseudo-thresholds. The presence of such higher-order 
poles can lead to large round-off errors, due to incomplete cancellations in computer arithmetic 
with finite precision. Fortunately, this is avoided in most cases by applying the partial fractions 
technique to rewrite the coefficients in the derivatives with respect to s, so that all poles in s are 
at most simple poles. This can always be done for the B, S, and T functions. 

For the U functions, double poles in the coefficient functions for the derivative with respect to s 
remain only if the second argument vanishes. Here, we take advantage of the facts that U(x, 0, y, z) 
does not enter into the differential equations that govern the other basis functions, and it can 
always be written algebraically in terms of them (see eqs. (A. 15), (A. 17) and (A. 18) of ref. [S]). 
Therefore, when the second argument of a U function vanishes, the result obtained for it from the 
Runge-Kutta running is irrelevant; it is simply replaced by the algebraic result before it is returned 
by the program, eliminating the roundoff error problem. 

In most cases for which double poles in the coefficient functions of the derivative of the master 
integral cannot be eliminated, it can be evaluated in terms of polylogarithms, so the Runge-Kutta 
technique is not needed anyway. A special case in which this does not occur is M(x, y, z, u, v) for 
v = {y/x + yjy) 2 ; then d(sM)/ds in ref. [J| has coefficients with double poles at sq = [\/x(u — y) + 
y/y(z — x)]/y/v. However, here we can use the identity: 

= yfx[U (z, x, y, v) - T(x, u, v)\ + y/y [U(u, y, x, v) - T(y, z, v)} 

+V^[T(v, u, x) + T(v, y, z)-l} + [A(v)fyfi - A(x)/yfx - s/y]B(y, u) 
+[A(v)/y/^-A(y)/y/y-y/x]B(x,z), (3.19) 

valid in general for v = (y/x + y/y) 2 - When the right side of equation (|3.19p is multiplied by 
[\fz(u — y) + y/u(z — x)]/v(s — so) 2 and added to the expression for d{sM)/ds from ref. [4], the 
result generically has no double poles in s. That is the form used by the program in this special 
case (and others related by permutations of the masses). 

A non-generic sub-case of the preceding, for which double poles in the coefficient functions of 
d{sM)/ds are not eliminated, is M(0, y, z, u, y) with y = (i/i± \fu) 2 . Because the coefficients 
involve double poles at s = z, there is some loss of accuracy at, and very near, that threshold. 
Fortunately, there is no good reason why a relation like y = {yfz ± \/u) 2 should hold exactly in 
a quantum field theory, unless a symmetry makes y or z or u equal to 0, and in each of these 
cases the master integral and all of its subordinates are given in terms of polylogarithms. More 
generally, we have checked that the coefficients in the derivatives as implemented in TSIL always 
have only simple poles, except in "unnatural" cases of this type (where no symmetry of a quantum 
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field theory can cause the necessary coincidence), or when the integrals are already analytically 
computed. 

The measures detailed above are generally sufficient to give good numerical accuracy near 
thresholds and pseudo-thresholds (except in the unnatural coincidence case just mentioned) without 
need for interpolation techniques or special expansions. 

4 Description of the program 

TSIL is a library of functions written in C, which can be (statically) linked from C, C++ or Fortran 
code[j] In addition to the main evaluation functions, it contains a variety of routines for I/O and 
other utilities. A complete list of functions in the user API is given in Appendix C. 

The principal data object in the code is a C struct with type name TSIL_DATA that contains 
values of the parameters x, y, z, u, v and Q 2 as well as the 15 basis functions of types B, S, T, U, M. 
Each integral function is itself a struct containing its value, arguments, and various unchanging 
coefficients used in computing its derivative. These coefficients are functions of x, y, z, u, v and are 
computed when the parameter values are set. In addition, each basis function contains a set of 
pointers to the other functions needed in evaluating its derivative. Also contained in the data struct 
are values of the integrals T, V, and S, T, U, and V. Definitions of all datatypes are contained in 
the header file tsil.h, which must be included in all application programs. 

In any program that calls TSIL functions requiring Runge-Kutta evaluation, at least one of 
these high-level data objects must be declared: 

TSILJDATA foo; 

(More than one such object, and arrays of such objects, are allowed. See subsection 15.31 for an 
example.) Users can of course access the items in the struct directly, though it is recommended 
that the provided user interface routines be used. These allow one to extract values of individual 
functions (or all of them), set the values of the external parameters, and so on. 

The size of the basic datatypes used for floating point values is controlled by the user (at compile 
time) with the choice of compiler flag -DTSIL_SIZE_LONG or -DTSIL_SIZE_DOUBLE. Then the type 
TSIL_REAL is accordingly synonymous with long double or double, while the type TSIL_COMPLEX 
is synonymous with long double complex or double complex. The recommended default size is 
long double on systems where it is available. For most physics applications (taking into account 
the natural suppression of two-loop effects), double should easily give sufficient accuracy. However, 
the use of long double provides a nice safety margin, and execution times are typically short (of 
order tenths or hundredths of a second on a modern workstation for generic inputs) in any case. 
Generally, long double data (typically with 63 or more bits of relative precision) gives results with 
relative accuracies better than 10 -10 for generic cases, but sometimes somewhat worse in cases with 
large mass hierarchies, and in some particularly difficult cases significantly worse. [The function 

wrapper routine is included that provides the interface to Fortran; see section [5741 below. 



13 



V(x,y,z,u) for very small but non-zero y can be particularly sensitive to roundoff errors, since 
the individual terms in its evaluation are proportional to 1/y and yet it is only logarithmically 
divergent as y — > 0.] The user should consider modifying the default parameters of the program if 
significant sensitivity to parameters is expected (or observed), or if speed is an overriding concern. 

In a typical application the user will initialize the data object and set values for the external 
parameters x,y, z,u,v,Q 2 by calling the function TSIL_SetParameters. (Parameter values may 
be changed in a data object with subsequent calls to this function.) Then the basis integrals 
are evaluated at any desired value of s by calling the master evaluation function TSIL jwaluate. 
Additional calls to TSIL_Evaluate can be used to re-compute the basis functions for other values 
of s. 

TSIL_Evaluate first decides whether the case at hand is known analytically; if so, the basis func- 
tions are computed directly. If not, numerical integration is required. In this case TSIL_Evaluate 
begins by rescaling all dimensionful quantities by the largest of x,y, z,u,v, |s|, so that all parame- 
ters are rendered dimensionless. It then locates all thresholds and pseudo-thresholds and decides 
whether special handling is needed, that is, if one of these special points is near s = or the 
final value of s. The nearness criterion is controlled by a constant THRESH_CUT0FF, defined in 
tsil_params .h. If s = is a threshold (or there is a threshold very near s = 0), evaluation pro- 
ceeds as described above by making the change of integration variable (|3.7p . If the final value of s 
is at or near a threshold, the change of variable (|3.8p is enabled for the final leg of the integration 
contour. 

Next, TSIL_Evaluate checks to see if the final value of s is smaller than the smallest non-zero 
threshold or pseudo-threshold; if so, then integration proceeds directly along the real s axis. If not, 
the generic displaced integration contour (fig. [2]) is used. In cases where the final destination s is 
near a threshold or pseudo-threshold, the 5-stage Runge-Kutta routine described by eqs. (|3.15p - 
(|3.18p is used for the very last step of the integration, to avoid evaluation of any derivatives at the 
endpoint. 

After the Runge-Kutta integration, the program checks to see if this was a case with uncanceled 
double pole terms in the derivatives of the U functions, due to the second argument vanishing. If 
so, these values are corrected. (Recall that in such cases the incorrect values obtained by the 
Runge-Kutta integration have no effect on other basis functions.) The program also replaces any 
subordinate integrals (B, S, T, U) that can be computed in terms of polylogarithms by their 
analytical values. Finally, the program computes the values of all T and V functions, and the 
"bold" variants of all functions defined in section 2, using eqs. (|2.35p - (|2,39p . 

The function TSIL_Evaluate returns 1 (TRUE) for successful execution or (FALSE) for error 
execution. A warning message is printed if the external parameters correspond to the unnatural 
threshold case discussed at the end of section 3. The data object further contains a status parameter, 
accessible via the function TSIL_GetStatus, which indicates how the master integral evaluation was 
performed: either analytic, numerical integration along real axis, or numerical integration along 
the contour of fig. [2 
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The standard output function is TSIL_PrintData, which prints all function values on stdout. 
An alternate format, designed so that captured output can serve as valid input files for Mathematica, 
is given by TSIL_PrintDataM. Additional utilities allow the user to extract individual basis functions 
or sets of functions to arrays. Note that warning and error messages appear on stderr so they 
may be redirected by the shell and examined separately. 

Along with the size of intrinsic datatypes, the parameters associated with the numerical in- 
tegration adaptive step-size control exert the main influences on execution speed and accuracy^] 
They are realized as members of the TSILJDATA struct, with names: 

• precisionGoal: This is 5p in eq. (|3.14|) . (We use a safety factor S = 0.9.) If the maximum 
estimated error for any dependent variable exceeds 5p multiplied by the increment of that 
variable for that step, and also exceeds the relative precision of the computer arithmetic times 
the absolute value of that variable, then the step is retried with a smaller step size, unless the 
step size would become smaller than specified below. Also, after a successful step, the size for 
the next step is chosen according to eq. (|3,14|) . unless it would exceed the amount specified 
below. (Defaults: 10~ 12 for long double, 5 x 1CT 11 for double.) 

• nStepsStart: For each leg of the contour, the initial step size is chosen so that there would 
be this many steps if the step size did not change. (Default: 500) 

• nStepsMin: The maximum allowed step size on a leg of the contour with dimensionless 
(rescaled) independent variable length L is given by L/nStepsMin. (Default: 100) 

• nStepsMaxCon, nStepsMaxVar: The minimum allowed step size on a leg of the contour with 
dimensionless independent variable length L is given by 

L/(nStepsMaxCon + L*nStepsMaxVar) . (Defaults: 10000, 10000) 

The step size is not allowed to increase by more than a factor of 1.5 or decrease by more than a factor 
of 2 after each step or attempted step. Note that by setting precisionGoal to 0, one can arrange 
that the total number of steps on each leg tends to nStepsMaxCon + L*nStepsMaxVar. If instead 
one sets precisionGoal to a very large number, the number of steps will tend to nStepsMin. The 
default values have been found to give good results for a wide variety of different choices of input 
parameters, for the integration variables used in the program. 

In addition to the evaluation for generic parameters described above, TSIL provides functions 
for direct analytical evaluation of the vacuum integrals A{x) and I(x,y,z), the one-loop inte- 
gral B(x,y), as well as B(x',y), dB(x,y)/ds, A e (x), B e (x,y), I(x',y,z), I(x",y,z), I(x',y',z), 
I{x"' ,y, z), and all S, T, T, U, V, M functions for which results in terms of polylogarithms are 
available in the literature. 

*These are always set, by TSIL_SetParameters, to be equal to the values specified at compile time in the file 
tsil_params .h. However, to deal with exceptional situations, they can optionally be reset at run time with the 
function TSIL_ResetStepSizeParams, after calling TSIL_SetParameters and before calling TSIL_Evaluate. 
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5 How to use the program 

5.1 Building TSIL 

Complete instructions for building the library are provided with the distribution. Typically the user 
edits the Makefile to choose the desired data size and set appropriate optimization flags for the 
compiler. The command make then produces the static archive libtsil.a, which may be placed 
in any convenient directory <dir>. The user then typically links to this library by passing the flags 

-lm -L<dir> -ltsil -DTSIL_SIZE_<size> 

to the linker, where the same flag -DTSIL_SIZE_<size> was used in the Makefile when compiling 
libtsil . a. The header file tsil .h must be included in any source file that makes use of the TSIL 
routines or data structures. 

The command make also produces an executable tsil, which takes as command-line arguments 
x, y, z, u, v, s, Q 2 and prints out the values of all integral functions together with timing and other in- 
formation. Also included with the package is a test program testprog . c (with executable tsil . tst 
produced by make) and a set of 320 files containing comprehensive test data. These include cases 
representing all known analytic results as well as cases requiring integration that have thresholds 
and pseudo-thresholds at s = and at the final s. Users should run this test suite after building the 
library to insure that accurate results are obtained. The test program uses pass/fail/warn criteria 
that are controlled by macros TSIUASS and TSIL_WARN in tsil_testparams .h. The first sets the 
maximum relative error allowed for the test to pass; the second sets a lower error threshold below 
which the test is deemed to fail. A relative error between these two values results in a warning. 

5.2 Essential functionality 

In the simplest application of TSIL, the parameters x, y,z,u,v > and Q 2 > will be set using 
TSIL_SetParameters, the integrals for real s evaluated using TSIL_Evaluate, and the results ex- 
tracted by the calling program with the command TSIL_GetFunction. The code for this might 
look like: 

TSIL_SetParameters (&foo, x, y, z, u, v, qq) ; 
TSIL_Evaluate (&foo, s) ; 

integrall = TSIL_GetFunction (&foo, <stringl>) ; 
integral2 = TSIL_GetFunction (&foo, <string2>) ; 

where foo has type TSIL_DATA, and x,y,z,u,v,s,qq all have type TSIL_REAL, and integrall, 
integral2, . . .have type TSIL_CDMPLEX, and <stringl>, <string2>, . . . can each be one of 

"M", "Uzxyv", "Uuyxv", "Uxzuv" , "Uyuzv" , "Tvyz", 

"Tuxv", "Tyzv", "Txuv" , "Tzyv" , "Tvxu" , "Svyz", "Suxv", 

according to which of the integrals in eq. (|2.40p is desired. The strings can also be one of 
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Vzxyv", "Vuyxv", "Vxzuv" , "Vyuzv" , "Bxz" , "Byu", 



to access the functions V and the one- loop B functions. Examples are given in subsection 15.31 The 
functions S, T, U and V can be accessed in a similar way, e.g.: 

integral3 = TSIL_GetBoldFunction (fefoo, <string3>, n) ; 

would return the coefficient of l/e n (for n = 0, 1, or 2) in the bold-faced function corresponding to 
an appropriate <string3> from the list above. 

All integrals that are analytically known in terms of polylogarithms can also be evaluated 
directly, without TSIL_SetParameters or TSIL_Evaluate or TSIL_GetFunction. For example, 

TSIL_Manalytic (x,y,z,u,v,s,&result) ; 

will return the int value 1 and set the variable result equal to M(x, y, z, u, v) for the appropriate 
s if it is analytically available, and otherwise will return 0. Here x,y,z,u,v are of type TSIL_REAL, 
and s, result must be of type TSIL_C0MPLEX. The functions TSIL_Sanalytic, TSIL_Tanalytic, 
TSIL_Tbaranalytic, TSILJJanalytic, and TSIL_Vanalytic have analogous behavior, except that 
they carry an additional argument qq of type TSIL_REAL for the renormalization scale squared Q 2 . 
For example, 

TSILJJanalyt ic (x , y , z , u , s , qq, feresult ) 

will return the int value 1 and set the variable result equal to U(x,y,z,u) for the appropriate 
s and Q 2 , if it is analytically available, and otherwise will return 0. The other analytic functions 
assign without pointers, for example 

result = TSIL_Bp (x,y,s,qq); 

will set result equal to B(x',y) computed analytically for the appropriate s and Q 2 . Some exam- 
ples are given in the next subsection, and a complete list of the TSIL user API is given in Appendix 
[6] and the program documentation and header files. 

In some applications, it could be that rather than a master integral M and all of its subordinates, 
one may only need integral functions corresponding to S, T, U or only S, T. Those cases can be 
evaluated simply by choosing any convenient numbers for the irrelevant squared-mass arguments. 
However, this is clearly not optimally efficient. In version 1.1 of TSIL, we have added the capability 
to only compute the integrals in the S, T, U system, or only those in the S, T system. This has been 
accomplished by adding functions TSIL_SetParametersSTU and TSIL_SetParametersST, each of 
which sets appropriate subsets of the squared mass parameters. A subsequent call of the function 
TSIL_Evaluate will only evaluate the relevant subset of integral functions. So, the user code might 
include for example: 

TSIL_SetParametersSTU (fefoo, x, z, u, v, qq) ; 
TSIL_Evaluate (fefoo, s) ; 
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integrall = TSIL_GetFunction (&foo, <stringl>) ; 
integral2 = TSIL_GetFunction (&foo, <string2>) ; 



where <stringl>, <string2>, . . . can each be one of 

"Uxzuv", "Tuxv", "Txuv", "Tvxu" , "Suxv", "Bxz" . 

Or, if U(x, z, u, v) is not needed, then one can use TSIL_SetParametersST (&foo, x, u, v, qq) ; 
instead. For generic cases, the S, T, U evaluation is a factor of 4 or 5 faster than full evaluation, and 
the S, T case gains a further 20% in evaluation speed. Note that the choice of which integral func- 
tions are evaluated by TSIL_Evaluate is controlled by the most recent call of TSIL_SetParameters 
or TSIL_SetParametersSTU or TSIL_SetParametersST for the relevant data struct. Also, if S, T, U 
or S, T evaluation is used, then TSIL_GetData and TSIL_GetBoldData will generate an error mes- 
sage; only TSIL_GetFunction and TSIL_GetBoldFunction should be used to extract the results in 
those cases. Note that in subset evaluation cases where there is only a single function of a given 
type (U, V, S, or B), the specification string may be abbreviated to only the first character, e.g. 
"U" in place of "Uxzuv" above. 

5.3 Sample applications 

As a sample application of TSIL, let us calculate the two-loop self energy and pole squared mass 
for a single scalar field with interaction Lagrangian 

£ = -ImV-|^ 3 -^ 4 . (5.1) 
Here m 2 , g and A are the tree- level renormalized parameters. The self-energy 

is a function of s = —p 2 , with p^ the external momentum, as well as the parameters m 2 , g and A. 
Note that the metric is either of signature ( — or Euclidean. Furthermore, s must be taken 
to be real with an infinitesimal positive imaginary part to properly resolve the branch cuts. 
The pole squared mass 

s = M 2 - iTM (5.3) 

is defined as the position of the pole, with non-positive imaginary part, in the propagator obtained 
from the perturbative Taylor expansion of the self-energy about the tree-level squared mass[j] (Note 
that in the present case the width T = identically.) This leads to the following iterative scheme 
for computing the two-loop pole squared mass. First, the one-loop approximation s^ 1 ' to the pole 
squared mass is obtained as 

S W = m 2 + — "-^nW (m 2 ) . (5.4) 
16ir z 



^In a theory with gauge fields the self-energy is gauge-dependent, but the pole squared mass defined in this way 
is properly gauge invariant. 
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Then, defining 

n 



16vr 2 



n«(m 2 ) + ( S « - m 2 )n«'(m 2 )] + ^J-^n( 2 )(m 2 ) , (5.5) 



in which the prime indicates a derivative with respect to s, we obtain the two-loop approximation 
to the pole squared mass as 

s ( 2 )=m 2 + n. (5.6) 
For the scalar theory described by eq. (|5.ip we find, including the MS counterterms: 

nW(i) = ^AA^-i^fox) (5.7) 

n( 2 )(s) = — — g 4 M(x, x, x, x, x) — — g 4 V(x, x, x, x) + g s U(x, x, x, x) 

--X 2 S(x,x,x) + \xg 2 B(x,x)B(x,x) + \x 2 A(x) [A(x)/x + l] (5.8) 
6 4 4 

-^Xg 2 A(x)B(x',x) - -Xg 2 I(x',x,x) , 



where x = m 2 , and 



rt 1 \s) = - 1 -g 2 ^B{x,x). 



(5.J 



Note that both il^ 1 )' and n( 2 ) have (1 - s/4x) _1 / 2 singularities at the threshold s = 4x, though this 
is not manifest in the above formulas. The TSIL code handles such true singularities by returning 
a value that is interpreted and displayed as the string Complexlnf inity. 

Below is a sample code that uses TSIL to calculate the pole squared mass for parameter values 
m 2 , g, X and Q 2 obtained as command-line inputs in that order. It first computes the required 
basis functions at s = m 2 , then assembles IlW(m 2 ), n( 2 )(m 2 ) and II"'(m z ), and finally outputs 
the pole squared mass: 

/* === scalarpole.c === 
* 

* Command-line arguments are: 

* scalar mass squared = x, 

* cubic coupling = g, 

* quartic coupling = lambda, 

* renormalization scale squared = qq. 
* 

* Run as, for example: ./spole 12 3 1 
*/ 

#include <stdio.h> 

#include "tsil.h" /* Required TSIL header file */ 

#define PI 4.0L*TSIL_ATAN(1 .OL) /* Uses arctan function defined in tsil.h */ 

int main (int argc, char *argv[]) 
{ 

TSIL_DATA result; /* Top-level TSIL data object */ 
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TSIL_REAL qq; /* Ensures correct basic type; see also TSIL_COMPLEX */ 

TSIL_REAL x, g, lambda; 

TSIL_CDMPLEX pil, pilprime, pi2, si, s2; 

TSIL_REAL factor = 1 . OL/ (16 . 0L*PI*PI) ; 

/* If incorrect number of args, print message on stderr and exit: */ 
if (argc != 5) 

TSIL_Error ("main" , "Expected 4 arguments: m~2, g, lambda, and Q~2", 1); 

/* Note cast to appropriate floating-point type for safety */ 
x = (TSIL_REAL) strtold(argv [1] , (char **) NULL); 

g = (TSIL_REAL) strtold(argv[2] , (char **) NULL); 

lambda = (TSIL_REAL) strtold(argv [3] , (char **) NULL); 
qq = (TSIL_REAL) strtold(argv [4] , (char **) NULL); 

/* All loop integrals have a common squared-mass argument x: */ 
TSIL_SetParameters (feresult, x, x, x, x, x, qq) ; 

/* For the pole mass calculation, evaluate two-loop integrals at s = x: */ 
TSIL_Evaluate (feresult, x) ; 

/* Assemble one- and two-loop mass squared results: */ 

pil = 0.5L*lambda*TSIL_A(x,qq) - . 5L*g*g*TSIL_B(x,x,x,qq) ; 

pilprime = -0.5L*g*g*TSIL_dBds(x, x, x, qq) ; 

pi2 = - 0.5L*g*g*g*g*TSIL_GetFunction (feresult, "M") 

- 0.5L*g*g*g*g*TSIL_GetFunction (feresult, "Vzxyv") 
+ g*g*g*TSIL_GetFunction (feresult, "Uzxyv") 

- (1.0L/6.0L)*lambda*lambda*TSIL_GetFunction(feresult, "Svyz") 

+ 0.25L*lambda*g*g*TSIL_P0W(TSIL_GetFunction(feresult, "Bxz"), 2) 
+ 0.25L*lambda*lambda*TSIL_A(x,qq)*(TSIL_A(x,qq)/x + 1.0L) 

- 0.5L*lambda*g*g*TSIL_A(x,qq)*TSIL_Bp(x, x, x, qq) 

- 0.25L*lambda*g*g*TSIL_I2p(x,x,x,qq) ; 

si = x + factor*pil; 

s2 = x + factor*pil + factor*f actor* (pi2 + pil*pilprime) ; 

printf ("Tree-level squared mass: °/ lf\n", (double) x) ; 
printf ("One-loop pole squared mass: 7 lf\n", (double) si); 
printf ("Two-loop pole squared mass: 7,lf\n", (double) s2) ; 

return 0; 

} 

Note the use of TSIL_A, TSIL_I2p, TSIL_B, TSIL_Bp, and TSIL_dBds to evaluate the functions 
A(x), I(x',x,x), B(x,x), B(x',x), and dB(x,x)/ds, respectively. [In the evaluation of pi2, we 
arbitrarily chose to use TSIL_GetFunction(&result , "Bxz") where TSIL_B could have been used.] 
The compile command for this program to produce the executable spole would typically be 

cc -o spole scalarpole.c -lm -L<dir> -ltsil -DTSIL_SIZE_<size> 

as discussed in subsection 15.11 
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Figure 3: Feynman diagram topologies for some two-loop corrections to the neutral Higgs scalar 
boson self-energies in supersymmetry. 

Other situations can be treated with appropriate generalizations. For example, in the Minimal 
Supersymmetric Standard Model correction to the neutral Higgs scalar boson self-energies |61| . 
there are graphs with the topology shown in fig. [3] involving the top quark t, the top squarks t n 
for n = 1, 2, and the gluino g. All of the necessary one-loop and two-loop basis integrals for these 
diagram topologies can be computed in one fell swoop with: 

TSILJDATA result [2] [2] ; /* Declare an array of TSIL_DATA structs */ 
TSIL_REAL mt2, mstop2 [2] , mgluino2; /* Squared masses of top, stops, gluino */ 
TSIL_REAL s; 
int i , j ; 

for (i=0; i<2; i++) { 
for (j=0; j<i; { 
TSIL_SetParameters (&(result [i] [j] ) ,mstop2 [i] ,mt2 ,mstop2 [j] ,mt2,mgluino2,qq) ; 
TSIL_Evaluate (& (result [i] [j]),s); 

} 

} 

The index of mstop2[] in the code is one less than the conventional top squark mass eigenstate 
label, because arrays start at index in C. Note that the evaluation for j = 1, i = would be 
redundant by symmetry with that for j = 0, i = 1, and so is not performed here. Note also that 
the subordinate integrals of topology S, T, and U are needed for the calculation of the self-energy 
and the pole mass, even though there are no Feynman diagrams with the corresponding topologies 
since there are no four-particle couplings involving fermions. The necessary two-loop integrals can 
then be extracted by, for example: 



value = 


TSIL_GetFunction(&(result [0] [0] ) , 


"M") ; 


for M(m| i , ml , 


mf, 


mg), 


value = 


TSIL_GetFunction(&(result [0] [0] ) , 


"Vuyxv") ; 


for V(mf ,mf , 


m|, 


m- g ), 


value = 


TSIL_GetFunction(&(result [1] [0] ) , 


"Vzxyv") ; 


for V(m~ , mf 


mf, 


mg), 


value = 


TSIL_GetFunction(&(result [1] [0] ) , 


"Vxzuv") ; 


for V(mk , mf 


2 


m-g), 


value = 


TSIL_GetFunction(&(result [0] [0] ) , 


"Txuv") ; 


for T(m~^ 


mf, 


m-g), 



etc. 



21 



5.4 Using TSIL with FORTRAN 

It is possible to call the TSIL functions from a Fortran program, and utilities for this are included 
with the package. Basic functionality is provided by a "wrapper" function tsilf evaluate, which is 
called as a subroutine from a Fortran program. This subroutine implements the most general TSIL 
calculation: it takes as arguments x,y, z,u,v,Q 2 , s and returns the values of all basis functions, 
including T, V, and "bold" functions. 

The results are returned to the calling program in a COMMON block, which corresponds to a 
special C struct used in tsilf evaluate. This COMMON block contains a number of pre-defined 
arrays that hold the various function values. Definitions of the COMMON block and subsidiary arrays 
are given in a header file that users include in their Fortran programs. In addition, this header file 
defines a set of integer variables that allow items in the COMMON block to be referred to by name. 

A Fortran program fragment that uses these utilities is shown below: 

PROGRAM ftest 
c Includes array and COMMON definitions: 

INCLUDE 'tsil_fort.inc' 

c (Code setting values of x,y,z,u,v,qq,s not shown) 

c Evaluate basis integrals: 

CALL tsilf evaluate (x,y,z,u, v,qq, s) 

c Print a representative value: 

PRINT *, U(xzuv) 

The provided wrapper code can serve as a model for users wishing to write their own interface 
routines with additional functionality. The TSIL documentation contains a detailed discussion of 
the issues that arise in using TSIL with Fortran. 

6 Outlook 

The TSIL library is intended to provide a convenient all-purpose solution to the problem of nu- 
merical evaluation of two-loop self-energy integrals in high-energy physics. The functions provided 
should work for arbitrary input parameters, without relying on special mass hierarchies or other 
simplifications that can arise in special cases. The library does take advantage of simplifications 
when they allow evaluation in terms of polylogarithms. 

The library is organized around the calculation of the basis integrals to which all other self- 
energy contributions can be reduced by known algorithms. It is therefore the responsibility of the 
user to perform the non-trivial calculations necessary to assemble the basis functions into physical 
observables. This involves the reduction of the Feynman diagrams to the basis integrals and the 
inclusion of counterterms, tasks that can always be automated using purely symbolic computer 
algebra manipulations or even performed by hand in favorable situations. We believe that keeping 
these separate from the problem of numerical evaluation is advantageous, and that the modular 
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nature of this approach will afford the flexibility to deal with the surprises that hopefully await us 
as we explore physics at the TeV scale. 

Appendix A: Comparison of conventions with other sources 

In this Appendix, we list the correspondence between our notation for the loop integrals and those 
found in several other references. In the following, the functions as defined in this paper and in 
refs. [I], [5] always appear on the left-hand side, and equivalent notations for them in other papers 
appear on the right. 

The definition of the master integral in ref. [16] is given by: 

M(x,y,z,u,v) = -I(s)/s, (A.l) 

with m\ , m\, m 2 , m% m\ = x, z, v, y, u, and the integral on the right-hand side defined in 4 dimen- 
sions. 

The notation used in refs. [3 EJ (HJ H3 Ej is: 

A(x) = -Aq(x) 

B(x,y) = B (s;x,y) 

I{x,y,z) = -T 135 (x,y,z) 

S{x,y,z) = -T 2 34(s;x,y,z) 

T(x,y,z) = T 22 34(s;x,x,y,z) 

U(x,y,z,u) = T 12 u{s;y,x,z,u) 

V(x,y,z,u) = -T 112 34(s;y,y,x,z,u) 

M(x,y,z,u,v) = -Ti234 5 (s;x,z,v,u,y). 

Ref. |36] used the same notation, with the exception that S(x,y, z) = —Ti23(s;x,y,z) there. 



The notation in ref. [T] is: 

A(x) = ial[ d) with m\ = x, (A. 10) 

B(x,y) = — iaG^ (s) with ml, m\ = x, y, (A. 11) 

I(x,y,z) = £i 2 J^q(0) with ml, m\, m\ = x, y, z, (A. 12) 

S(x,y,z) = a 2 j[f[(s) with m\ , m\, m\ = x, y, z, (A. 13) 

T(x,y,z) = -a 2 J^\(s) with m\, m\, m\ = x, y, z, (A. 14) 

XJ(x,y,z,u) = —a 2 Viil 1 (s) with m\ , m\, m\, m\ = u, x, z, y, (A. 15) 

\{x,y,z,u) = a 2 Viil 2 (s) with m\, m\, m\, m\ = u, x, z, y, (A. 16) 

M(x, y, z, u, v) = a 2 Fifl 11 (s) with m\ , m\, m\, m\, m\ = x, y, z, u, v , (A. 17) 

and the closely related notation of ref. [2] is 

A(x) = mTAl[d,s,{{l, y^}}], (A.18) 
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(A.2) 
(A.3) 
(A.4) 
(A.5) 
(A.6) 
(A.7) 
(A.8) 
(A.9) 



B(x,y) 
l[x,y,z) 
S(x,y,z) 
T(x,y,z) 
U(x,y,z,u) 
V(x,y,z,u) 
M(x, y, z, u, v) 



-iaTBl[d,s,{{l,^,{l,^y}}] 
a 2 TJl[d,0,{{l,^},{l,^y},{l,^I}}] 

a 2 TJl[d,s,{{l,^},{l,^y},{l,^}}] 
-a 2 TJl[d, s ,{{2, y/x~}, {1, yy}, {1, ^}}] 
-a 2 TVl[d, s , {{1, {1, V^}, {1. v^}, {1, ^y}}] 
a 2 TVl[d, s , {{1, y/u\, {1, V^}, {1, v^K {2, Vy}}] 
a 2 TFl[d, s ,{{1, v^}, {1, ^y}, {1, {1, V^}, {1, v^}}] 



(A.19) 
(A.20) 
(A.21) 
(A.22) 
(A.23) 
(A.24) 
(A.25) 



where a = (Attu 2 ) 2 d l 2 . 

The notation of [32j [531 ES [56j [57] (ref. [5l] used a slightly different normalization) is: 



A(x) = 


bT(d,x) 






B(x,y) = 


bS(d, x, y, - 






Kx,y,z) = 


b 2 V(d,x,y, 


*) 




S(x,y,z) = 


b 2 F (d,x,y 


z, 


s) 


T(x,y,z) = 


b 2 Fi(d,x,y 


z, 


-s) = b 2 F 2 (d,z,x,y 


XJ(x,y,z,u) = 


b 2 F^(d, x, u 


y, 


z,-s), 


where b = Afi 4 ~ d . 








The notation of [50] is: 








S(x,y,z) 


= c 2 S lu 




with m\ , 777,2, m 3 = 



(A.26) 
(A.27) 
(A.28) 
(A.29) 
(A.30) 
(A.31) 



U(x,y,z,u) 
M(x,y,z,u,v) 



c 2 S 121 



2 r<221 



c z S 



x,y,z, 

with m\ , m 2 ,m 2 , m\ = z, u, y, x, 

■i-u 2 2 2 2 2 
With 771]^, 7772, 7773, 777-4, m 5 = x ; z > u > 2/) u > 



(A.32) 

(A.33) 
(A.34) 



with c= (2vr) ,i - 4 . 

Appendix B: Analytic expressions for some special cases 

In this Appendix, we present some analytic formulas for some important two-loop basis integrals 
special cases, two of which do not seem to have appeared before in the literature, and others that 
are equivalent to known results. 

As a generalization of the integral M(0, 0, x, x, 0) computed in ref. [16] . we find that for x > y: 



M(0,0,x,y,0) 



3Li 3 (r x ) + 3Li 3 (r J/ ) - 3Li 3 (r y /r x ) - 3Li 3 (yr y /xr x ) 
+3Li 3 (7//x) - 3C(3) + [hx(r y ) - ln(r x ) + 2 \n{y/x)]L\ 2 {yr y /xr x ) 
+ [ln(r y ) - ln(ra.) - 2\n(y / x)]L\ 2 (y / x) + \hx(r y ) - 3ln(r x )]Li 2 (r x ) 
+ [ln(r x ) - 3\n(r y )]Li 2 (r y ) + S[hx(r y ) - \n{r x )]L\ 2 (r y /r x ) 
+ [ln(r x ) - ln(r y )][ln(r x ) - ln(r y ) + ln(x/y)] ln(l - y/x) 
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M(0,x,y,0,y)\ ( 



-ln 3 ^) + 2\n 2 (r x )\n(r y ) - \n(r x ) \n 2 (r y ) - ln(r x ) ln(r y ) ln(l - r x ) 
+ ln(r x )ln(r y )ln(y/x) -la(r x )ln(y/x)[ln(r x ) +ln(y/x)]/2 
-C(2)[31n(r as ) + ln(r 1 ,)]]/s (B.l) 

where r x = 1 — s/x and r y = 1 — s/y are each given an infinitesimal negative imaginary part. The 
result for x < y is obtained by interchanging the squared masses. 

Generalizing the result M(0,x,x,0,x)\ s=x found in ref. [16], we obtain threshold cases: 

M(0,y, y,0, x)\ s=x = [2Li 3 ([r - l]/[r + 1]) - 2Li 3 ([l - r]/[r + 1]) - 2Li 3 (r/[r + 1]) 

-2Li 3 (l - 1/r) - Li 3 (l/r 2 )/4 + 2[ln(2r) - ln(r + l)]Li 2 (r/[r + 1]) 
+ [41n(r - 1) - 21n(2r) - 21n(r + l)]Li 2 (l - 1/r) 
+2[ln(r - 1) - ln(r + l)][Li 2 ([l - r]/2) - Li 2 ([r - l]/2r) 
- ln(4r)Li 2 (l/r 2 )/2 + ln 2 (r)[31n(r - 1) - ln(r + 1) - 61n2]/2 
+ ln(r)[21n(r + 1) ln(2[r - 1]) - ln 2 (r - 1) - 6C(2)] - (4/3) ln 3 (r) 
-(2/3)ln 3 (r + 1) +ln21n 2 (r + 1) +C(3)/2 + 6C(2)ln(l +r)j/x, (B.2) 

2Li 3 ([v^ - + 1]) " 2Li 3 ([l - v^]/[v^ + 1]) 

+2Li 3 (l - r) - 2Li 3 ([r - l)/r) - 3C(3)/2 
+ [ln(r) - 21n(r - l)]Li 2 (l - r) - [ln 2 (r - 1) ln(r)]/2 



+ [ln 3 (r)]/3 - 2C(2) ln(r) + 6C(2) ln(l + ^F) /a 



(B.3) 



where r = y/x is given an infinitesimal negative imaginary part to get the correct branches. 
Eq. (|B.3h is equivalent, by use of the recurrence relations of [I], to results already obtained in 
[361 EH- 

We also note the following threshold and pseudo-threshold values (equivalent to results obtained 
in |38|: see also |39|): 



U(x,0,y,y)\ s=x 

where r = \/yJx, and 
U(y,0,y,x)\ s=x = 



11/2 - 31nx + lnxlny - (lm/) 2 /2 + (1 + y/x){((2) - Li 2 (l - x/y)] 
-4r{Li 2 ([l - r]/[l + r)) - Li 2 ([r - l]/[r + 1]) + 3C(2)/2} (B.4) 



11/2 - 21nx - \ny + (lny) 2 /2 - 2C(2)(1 + y/x) 

+(rna; - 1 + y/x[l - lay]) ln(l - x/y - is) + (1 + 2y/x)Li 2 (l - x/y). (B.5) 



Appendix C: The TSIL Application Programmer Interface 

This Appendix lists the functions in the TSIL package, and their basic functionality. Complete 
details may be found in the TSIL documentation and header files. 
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Basic evaluation functions: 
TSIL_SetParameters 

TSIL_SetParametersSTU 

TSIL_SetParametersST 

TSIUvaluate 

TSIL_GetStatus 

TSIL_GetData 

TSIL_GetBoldData 

TSIL_GetFunction 

TSIL_GetBoldFunction 



Sets parameters x, y, z, u, v, Q 2 and selects evaluation of all 
integral functions 

Sets parameters x, z, u, v, Q 2 and selects evaluation of S, T, U 

functions only (new in vl.l November 2006; see end of section 5.2) 

Sets parameters x, u, v, Q 2 and selects evaluation of S, T 

functions only (new in vl.l November 2006; see end of section 5.2) 

Evaluate integral functions for a specified s 

Return current evaluation status 

Extract a set of integral function values to an array 

Extract a set of bold integral function values to an array 

Return a single specified integral function value 

Return a single specified bold integral function value 



I/O related functions: 
TSIL_PrintStatus 
TSIL_PrintData 
TSIL_WriteData 
TSIL_PrintDataM 
TSIL_WriteDataM 
TSIL_cprintf 
TSIL_cprintfM 
TSIL_Error 
TSIL_Warn 
TSIL_PrintInf o 



Print evaluation status to stdout 

Print all integral function values to stdout 

Write all integral function values to a file 

As TSIL_PrintData, but format is valid Mathematica input 

As TSIL_WriteData, but format is valid Mathematica input 

Generic printing of values of TSIL_Complex type 

As TSIL_cprintf , but in Mathematica input form 

Print a message to stderr and exit 

Print a warning message to stderr 

Write general information on stdout 



Utilities: 

TSIL_ResetStepSizeParams 
TSIL_IsInf inite 
TSIL_DataSize 
TSIL_PrintVersion 



Sets new parameters used for Runge-Kutta step size control 
Tests whether argument of type TSIL_Complex is finite 
Returns size of intrinsic floating point data used 
Prints version number of TSIL 



Analytic cases: 

TSIL_Dilog 

TSIL_Trilog 

TSIL_A 

TSILJVeps 

TSIL_B 

TSIL_Bp 

TSIL_dBds 

TSIL_Beps 



Dilogarithm function of complex argument 1^2(2) 

Trilogarithm function of complex argument Li3(z) 

One-loop vacuum integral A(x) 

One-loop vacuum integral A e (x) 

One-loop self-energy integral B(x,y) 

B(x',y) 

dB(x, y)/ds 

One-loop self-energy integral B e (x,y) 
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TSIL_I2 

TSIL_I2p 

TSIL_I2p2 

TSIL_I2pp 

TSIL_I2p3 

TSIL_S analytic 

TSIL_Tanalytic 

TS IL_Tbar analy t i c 

TSILJJanalytic 

TSIL_Vanalytic 

TSILJVI analytic 



Two-loop vacuum inte. 
I(x',y,z) 
I(x",y,z) 
I{x',y',z) 
I(x'",y,z) 

Analytic evaluation of 
Analytic evaluation of 
Analytic evaluation of 
Analytic evaluation of 
Analytic evaluation of 
Analytic evaluation of 



gral I(x, y, z) 



S if available 
T if available 
T if available 
U if available 
V if available 
M if available 



Fortran interface function: 



tsilf evaluate_ Wrapper for TSIL_Evaluate, callable from Fortran 
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